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Abstract 

Two new versions of the k — u two-equation turbulence 
model will be presented. The new Baseline (BSL) model 
is designed to give results similar to those of the original 
k — u model of Wilcox, but without its strong dependency 
on arbitrary freestream values. The BSL model is identical 
to the Wilcox model in the inner 50% of the boundary-layer 
but changes gradually to the standard k - e model (in a 
k - u formulation) towards the boundary-layer edge. The 
new model is also virtually identical to the k - e model for 
free shear layers. The second version of the model is called 
Shear-Stress Transport (SST) model. It is a variation of the 
BSL model with the additional ability to account for the 
transport of the principal turbulent shear stress in adverse 
pressure gradient boundary-layers. The model is based on 
Bradshaw’s assumption that the principal shear-stress is pro- 
portional to the turbulent kinetic energy, which is introduced 
into the definition of the eddy-viscositv. Both models are 
tested for a large number of different flowfields. The results 
of the BSL model are similar to those of the original k - u 
model, but without the undesirable freestream dependency. 
The predictions of the SST model are also independent of 
the freestream values but show better agreement with exper- 
imental data for adverse pressure gradient boundary -layer 
flows. 

Introduction 

The main field of application of Navier-Stokes meth- 
ods in aerodynamics will be for complex turbulent flows 
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that cannot be treated by inviscid, or viscous-inviscid inter- 
action schemes. Examples are massively separated flows, 
flows involving multiple length-scales, flows with three- 
dimensional separation and complex unsteady flows. In 
these flows, the application of algebraic turbulence mod- 
els, like the Cebeci-Smith [1], the Baldwin-Lomax [2] or 
the Johnson-King [3] model, becomes very complicated and 
often ambiguous, mainly because of the difficulty to define 
an algebraic length-scale. It is obvious that the improve- 
ment of numerical methods must be accompanied by the 
development of more general turbulence models and their 
implementation into Navier-Stokes codes. 

In addition to being independent of the specification of 
an algebraic length-scale, there is a long list of characteristics 
a good turbulence model wouM have to satisfy. Obviously, 
the model should be "sufficiently" accurate for the intended 
type of applications. Furthermore, and almost as important, 
the model must be numerically robust and should not con- 
sume excessive amounts of computation time (compared to 
the mean-flow solver). Another important demand is that the 
results should not have a strong dependency on ambiguous 
quantities, especially on the specified freestream values. 

The most popular non-algebraic turbulence models are 
two-equation eddy-viscosity models. These models solve 
two transport equations, generally one for the turbulent ki- 
netic energy and another one related to the turbulent length- 
(or time-) scale. Among the two-equation models, the k - e 
model is by far the most widely used today. The first low 
Reynolds number k -,e model was developed by Jones and 
Launder [4] and has subsequently been modified by many 
authors. 

The k - e model has been very successful in a large 
variety of different flow situations, but it also has a num- 
ber of well known shortcomings. From the standpoint of 
aerodynamics, the most disturbing is the lack of sensitivity 
to adverse pressure-gradients. Under those conditions, the 
model significantly overpredicts the shear-stress levels and 
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thereby delays (or completely prevents) separation [5]. Rodi 
[6] attributes this shortcoming to the overprediction of the 
turbulent length-scale in the near wall region and has shown 
that a correction proposed by Hanjalic and Launder improves 
the predictions considerably. However, the correction is not 
coordinate-invariant and can therefore not be applied gener- 
ally. An alternative way of improving the results has been 
proposed by Chen and Patel [7] and by Rodi [8], They re- 
place the e-equation in the near wall region by a relation 
that specifies the length-scale analytically. This also reduces 
some of the stiffness problems associated with the solution 
of the model. Although the procedure is coordinate indepen- 
dent, it has only been applied to relatively simple geometries, 
where the change between the algebraic relation and the e - 
equation could be performed along a pre-selected gridline. 
Clearly this cannot be done in flows around complex geome- 
tries. Furthermore, the switch has to be performed in the 
logarithmic part (the algebraic length-scale is not known in 
the wake region), so that the original k — t model is still 
being used over most of the boundary layer. 

Another shortcoming of the k - e model is associated 
with the numerical stiffness of the equations when integrated 
through the viscous sublayer. This problem clearly depends 
on the specific version of the k - e model selected, but there 
are some general aspects to it. All low Reynolds number k - e 
models employ damping functions in one form or another in 
the sublayer. These are generally nonlinear functions of di- 

l2 

mensionless groups of the dependent variables like Rf = . 

The behavior of these functions cannot easily be controlled 
by conventional linearization techniques and can therefore 
interfere with the convergence properties of the scheme. A 
second problem is that e does not go to zero at a no-slip sur- 
face. That in turn leaves two alternatives. One is to employ 

a nonlinear boundary condition on e (e — /(~^))> or t0 
add additional terms to the c - equation [4] that allow the use 
of a homogeneous boundary condition. Both approaches in- 
troduce additional nonlinearities that can upset a numerical 
procedure. 

There is a significant number of alternative models 
[9, 10, 1 1] that have been developed to overcome the short- 
comings of the k - e model. One of the most successful, 
with respect to both, accuracy and robustness, is the k — u 
model of Wilcox [9]. It solves one equation for the turbulent 
kinetic energy k and a second equation for the specific turbu- 
lent dissipation rate (or turbulence frequency) w. The model 
performs significantly better under adverse pressure-gradient 
conditions than the k - e model although it is the author’s 
experience that an even higher sensitivity to strong adverse 
pressure-gradients would be desirable [12]. Another strong- 
point of the model is the simplicity of its formulation in the 
viscous sublayer. The model does not employ damping func- 
tions and has straightforward Dirichlet boundary conditions. 
This leads to significant advantages in numerical stability. 


However, the k - u model also has a major shortcom- 
ing. It has been reported recently that the results of the 
model depend strongly on the freestream values, w y, that 
are specified outside the shear-layer. In [13] this problem 
has been investigated in detail, and it has been shown that 
the magnitude of the eddy-viscosity can be changed by more 
than 100% just by using different values for This is 
clearly unacceptable and corrections are necessary to ensure 
unambiguous solutions. 

In this paper, two new turbulence models will be pre- 
sented. First, a new baseline (BSL) k - w model will be 
described. It is identical to the k - w model of Wilcox [9] 
for the inner region of a boundary layer ( up to approximately 
6 /2) and gradually changes to the standard k - e model in the 
outer wake region. In order to be able to perform the com- 
putations with one set of equations, the k - e model was first 
transformed into a k - w formulation. The blending between 
the two regions is performed by a blending function Fy that 
gradually changes from one to zero in the desired region. No 
a priori knowledge of the flowfield is necessary to perform 
the blending. The function also ensures that the k - e for- 
mulation is selected for free shear-layers. The performance 
of the new (BSL) model is very similar to that of the Wilcox 
k — w model for adverse pressure gradient boundary-layer 
flows (and therefore significantly better than that of the k - e 
model), but without the undesirable freestream dependency. 
For free shear layers the new model is basically identical 
to the k — ( model, which predicts spreading rates more 
accurately than the k — w model. 

Although the original - and the new BSL k — w mode! 
perform better in adverse pressure gradient flows than the 
k - e model, they still underpredict the amount of separa- 
tion for severe adverse pressure gradient flows [12]. In an 
attempt to improve matters, the eddy-vicosity formulation of 
the BSL model will be modified to account for the transport 
effects of the principal turbulent shear-stress. The motiva- 
tion for this modification comes from the Johnson-King (JK) 
model [3] which has proven to be highly successful for ad- 
verse pressure gradient flows. The JK-model is based on 
the assumption that the turbulent shear-stress is proportional 
to the turbulent kinetic energy in the logarithmic and wake 
regions of a turbulent boundary layer. Johnson and King 
solve an equation for the maximum turbulent shear-stress 
at each downstream station and limit the eddy-viscosity in 
order to satisfy this proportionality. In the framework of 
two-equation models the turbulent kinetic energy is already 
known and it is therefore only necessary to limit the eddy- 
viscosity to account for the same effect. The resulting model 
will be called shear-stress transport (SST) model. First pre- 
dictions based on this assumption have already been reported 
in [12], 

The BSL and the SST models are not significantly more 
complicated than the original k-u> model and consume only 
little more computing time. Because the two models are 
virtually identical to the original k — u model in the near wall 
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region, the modifications developed by Wilcox [9] for rough 
walls and for surface mass injection can be applied without 
changes. Furthermore, the models have shown the same 
numerical robustness as the original model for all the flows 
computed so far. The present paper is based on the work 
presented in [14], However the equations have been changed 
somewhat in order to optimize the model performance for 
transonic flows and through transition. 

The Turbulence Model 
The new Baseline (BSL) Model 
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Tlie idea behind the BSL model is to retain the robust 
and accurate formulation of the Wilcox k - u model in the 
near wall region, and to take advantage of the freestream 
independence of the k - e model in the outer part of the 
boundary-layer. In order to achieve this, the £ - r model 
is transformed into a k - u formulation. The difference 
between this formulation and the original k - w model is that 
tin additional cross-diffusion term appears in the w-equation 
and that the modeling constants are different. The original 
model is then multiplied by a function F\ and the transformed 
model by a function ( 1 - Fj ) and both are added together. 
The function Fj will be designed to be one in the near wall 
region (activating the original model) and zero away from 
the surface. The blending will take place in the wake region 
of the boundary-layer. The left hand side of the following 
equations is the Lagrangian derivative: Dj Dt := d/dt + 
ujd/dx j. 

Original k — w model: 


Let o\ represent any constant in the original model (cr k \ , ...), 
<i>2 any Constantin the transformed k - e model (cr^, ...)and 
<p the corresponding constant of the new model (<x ;. . .), then 
the relation between them is: 

d = 0) 

the following two sets of constants will be used: 

Set 1 ( 0 j) (Wilcox): 

ff*i =0.5,<r wl =0.5, ft =0.0750, ( 8 ) 

P* = 0.09, k = 0.41, 7l = 0 l /0* ~ a jj\ K ~ i \f0* 

Set 2 (< fa ) (standard k - e): 

<r kl = 1 . 0 , tr u2 = 0-856, 0 2 = 0.0828, (9) 

0* = 0.09, k = 0.41, 72 = 02/0* ~~ ( 7 w 2 K ~/\/^*- 
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Transformed k - e model: 


Set 1 corresponds to the original k - u model and will be 
used in the near wall region exclusively. Set 2 corresponds 
to the transformation of the standard k - t model (q c = 
1 .44, ct £ = 1 .92) and its main area of application is for free 
shear-layers. 

The model has to be supplemented by the definition of 
the eddy-viscosity: 
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Now, equation (1) and equation (2) are multiplied by F x and 
equation (3) and equation (4) are multiplied by (1 — Fj) and 
the corresponding equations of each set are added together 
to give the new model: 
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The turbulent stress tensor r tJ = -pu-uj is then given by: 
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In order to complete the derivation of the model it is 
necessary to define the blending function F x . Starting from 
the surface, the function should be equal to one over a large 
portion of the boundary layer in order to preserve the de- 
sirable features of the k - w model, but go to zero at the 
boundary layer edge to ensure the freestream independence 
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of the k - f model. The function will be defined in terms of 
the variable: 


argi = min(mai :( 
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as follows: 

F\ = tanh(arg^) (13) 

where CD ^ is the cross-diffusion term of equation (4): 


equation for the turbulent shear-stress r that is based on 
Bradshaw’s assumption that the shear-stress in a boundary- 
layer is proportional to the turbulent kinetic energy, k, : 

t = pa\k (16) 

with a j being a constant. On the other hand, in two-equation 
models, the shear-stress is computed from: 


t - n t n 


(17) 




(14) 


with Q = For conventional two-equation models this 
relation can be rewritten to give: 


The first argument in equation (12) is the turbulent length 
scale Lf = '/k/( 0.09w)(= k 2 / 2 /(), divided by the shortest 
distance to the next surface, y. The ratio Lt/y is equal 
to 2.5 in the logarithmic region of the boundary-layer and 
goes to zero towards the boundary-layer edge. The second 
argument in equation (12) ensures that the function F\ does 
not go to zero in the viscous sublayer. The third argument 
is an additional safeguard against the “degenerate” solution 
of the original k — ui model with small freestream values 
for w [13]. Figure la shows the typical behavior of the 
function Fj for different velocity profiles in a strong adverse 
pressure gradient boundary layer (it also depicts the function 
F 2 explained later). Figure 1 also includes the underlying 
velocity profiles (same linestyle). The function is equal to 
one over about 50% of the boundary-layer and then gradually 
goes to zero. The behavior of the new BSL model will 
obviously lie somewhere between the original k-w and the 
k - t model. However, since most of the production of both, 
k and w. takes place in the inner 50% of the layer, it can be 
expected that the model performance will be closer to that 
of the k - u model, governing this area. Recall that the 
replacement of the e equation by an algebraic length-scale, 
as proposed by [7, 8] has to be performed in the logarithmic 
region so that the original k — e model still covers the largest 
part of the boundary layer and results will therefore be much 
closer to those of the k - e model. 


The Shear-Stress Transport (SST) Model 

One of the major differences between eddy-viscosity 
and full Reynolds-stress models, with respect to aerodynamic 
applications, is that the latter accounts for the important 
effect of the transport of the principal turbulent shear-stress 
r =: -pu'v 1 (obvious notation) by the inclusion of the term 


Dt _ dr 'dr 

Dt dt dxf. 
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The importance of this term has dearly been demonstrated 
by the success of the Johnson-King (JK) [3] model. Note 
that the main difference between the JK - model and the 
Cebeci-Smith model lies in the inclusion of this term in the 
former, leading to significantly improved results for adverse 
pressure gradient flows. The JK model features a transport 
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as noted in [12]. In adverse pressure gradient flows the ratio 
of production to dissipation can be significantly larger than 
one, as found from the experimental data of Driver [15], and 
therefore equation (18) leads to an overprediction of r. In 
order to satisfy equation (16) within the framework of an 
eddy-viscosity model, the eddy-viscosity would have to be 
redefined in the following way: 


The rational behind this modification can also be explained in 
the following way: In conventional two-equation models the 
turbulent shear-stress responds instantaneously to changes in 
the shear-strain rate Q. much likean algebraic eddy- viscosity 
model, whereas equation (19) guarantees that r does not 
change more rapidly than pa^k. Obviously, equation (19) 
is not desirable for the complete flowfield, since it leads to 
infinitely high eddy-viscosities at points where f2 goes to 
zero. Note however, that in adverse pressure gradient flows, 
production is larger than dissipation for the largest part of 
the layer (or Q > aj w). The following expression: 


max(ajw;n) 

guarantees therefore the selection of equation (19) for most 
of the adverse pressure gradient regions (wake region of the 
boundary layer), whereas the original equation (10) is used 
for the rest of the boundary layer. Figure 2 shows the relation 
of (— uV/ai k) versus \J Production/ Dissipation for the 
SST model (equation (20)), the conventional k — w (k - e) 
model (equation (10)), Bradshaw’s relation (equation (16)) 
and a relation proposed by Coakley [ 10] . Note that Coakley ’s 
relation contains the relations of Bradshaw (/? = 1) and that 
of the conventional two-equation models (/? = 0) as a subset, 
but not equation (20) (0 is a free parameter of Coakley’s 
model). 

In order to recover the original formulation of the eddy- 
viscosity for free shear-layers (where Bradshaws assump- 
tion, expressed in equation (16) does not necessarily hold) 
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the modification to the SST model has to be limited to wall 
bounded flows. This can be done in the same way as for the 
BSL model by applying a blending function F 2 . 


i / 

maxfapw; SlF 2 ) 

where F 2 is defined similarly to equation (13): 

. y/k 500l/ . 

ar 92 = mfiX ( 2 FTnQ — : “T - ) ( 22 ) 

0.09wy 

F 2 = tanh(arg2) (23) 

Fi is depicted in Fig. lb in the same way as Fj in Fig. la. 
Since the modification to the eddy-viscosity has its largest 
impact in the wake region of the boundary layer, it is im- 
perative that F 2 extends further out into the boundary-layer 
than Fi . (Note on the other hand that Fj has to fall off to 
zero well within the boundary-layer in order to prevent the 
freestream dependence of the k - u model). 

This modification to the eddy-viscosity is used in con- 
nection with the BSL model derived above. Flowever, in 
order to recover the correct behavior for a flat plate boundary 
layer, the diffusion constant had t0 be adjusted accord- 
ingly. The constants for the SST model are: 

Set 1 (SST - inner): 

<r kl = 0.85,<r wl =0.5, ft = 0.0750,0! =0.31, (24) 

6* = 0.09, k = 0.41, 71 = ft/,5* - <r ul K 2 /Jj3*. 

Set 2 remains unchanged. Furthermore, for general flows Q 
is taken to be the absolute value of the vorticity. Both models 
are given in their full form in Appendix A. 

Two new turbulence models have been introduced in 
this section. Both can be regarded as zonal models, since 
they utilize different models for different areas of the flow- 
field. However, in contrast to the classical zonal approach 
the present models do not require an a priori knowledge of 
the flowfield in order to define the zonal boundaries where 
the different models are to be used. The change between 
the different sub-models is achieved by "smart" functions 
that can distinguish between the different zones. Note that 
similar functions could also be designed for the k - e model. 
This makes it possible to design one set of constants for wall 
bounded flows and a second set for free shear-layers and 
to switch between the different sets in the same way as in 
equation (7). 

The BSL model retains the advantages of the Wilcox 
k — w model for boundary-layer applications, but avoids its 
undesirable freestream dependency. Furthermore it switches 
to the more accurate k - t model for free shear-layer ap- 
plications. In addition to this, the SST model modifies the 


definition of the eddy-viscosity for adverse pressure gradient 
boundary-layer flows in much the same way as the lohnson- 
King model does. From a computational point of view both 
models are not significantly more complex than the original 
k — u) model. 

Boundary Conditions 

At a no-slip wall, all turbulent quantities, except w are 
set to zero. As pointed out by Wilcox [9], w satisfies the 
following equation near the wall: 

6 // 

w — - — 7 as y -* 0. (25) 

ft 

Wilcox recommends to specify this analytical solution for 
the first few gridpoints away from the wall explicitly. The 
present author found it much easier and as accurate to im- 
plement the following boundary condition: 

6u 

uj = 10 x at y = 0 (26) 

ft (Ay) 2 

where Ay is the distance to the next point away from the wall. 
Equation (26) simulates the boundary condition 25 without 
the need of changing the solution at interior points. It should 
be noted that models based on the w-equation give accurate 
results if the near wall values of w are sufficiently large. Both, 
equations (25) and (26) satisfy this demand. The results are 
not sensitive to the factor (10) used in equation (26). 

At inflow boundaries, the turbulence quantities are spec- 
ified and at outflow boundaries a zero gradient is assumed. 

Two of the computed flowfields have a rotational sym- 
metry. In this case, the gradients of all turbulence quantities 
in the circumferential direction are set to zero. 

Numerical Method 

The mean flow equations are solved by the 1NS3D code 
of S. E. Rogers and D. Kwak [ 16] which is based on a pseudo- 
compressibility method. The convective terms are upwind 
differenced with a third-order flux-difference scheme. The 
viscous fluxes are differenced with second-order accurate 
central differences. The linear equations resulting from the 
first-order backward Euler time differencing are solved with 
an implicit line relaxation scheme. 

The turbulence equations have been solved with a num- 
ber of different schemes, from a third-order upwind differ- 
encing, to a second-order TVD (total variation diminishing) 
[17] to a first order upwind scheme. It was found in all 
computations that the solutions were virtually independent 
of the scheme used in the turbulence equations, although a 
change of accuracy in the mean flow solver had a large ef- 
fect. The reason for the insensitivity to the treatment of the 
convection terms in the turbulence model is, that unlike in 
the Navier-Stokes equations, they are not the leading order 
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terms. For this reason most of the computations are based on 
a first order upwind scheme. The turbulence equations were 
solved decoupled from the mean flow equations in a separate 
subroutine. One of the important aspects in the discretization 
of turbulence models is the implicit treatment of the source 
terms. The following approximate linearization gave good 
numerical properties: 


All computations have been performed on different 
grids, to ensure that the presented solutions are grid inde- 
pendent. The airfoil computations were performed on a 
standard grid kindly provided by S. E. Rogers [18]. 

Results 


(27) 

(P* + C u - D u ) sa - |C ^ 1 + 1Du ■ (28) 

OU) u 

where P, D. C are the production, the destruction and the 
additional cross diffusion terms respectively. The above 
expressions go to the left hand side of the algorithm with a 
change of sign and thereby increase the diagonal dominance 
of the method. The resulting numerical scheme has proven 
to be very robust and all of the computations with the k - u 
models could be run with an infinite time step. An exception 
is the backward facing step flow where the time step had 
to be reduced for all models tested so far by the author. 
Furthermore, the computations could be started with very 
crude initial conditions (like freestream values). 

Experience with two-equation turbulence models has 
shown that in regions of low values of u (f/fc), small dis- 
turbances in the shear strain rate can lead to erroneous 
spikes in the eddy-viscosity in the freestream or near the 
boundary layer edge. In order to understand the effect, 
the transport equation for the eddy-viscosity was derived 
from the k - u model for incompressible flows (note that 
Dvt/Dt = 1 /uDk/Dt - k/u l Du/Dt): 




^t/du ,• du,du x 


Dt 


' dxj dx{ dxj 


Flat Plate Boundary Layer 

In order to show the motivation for the derivation of the 
BSL model, flat plate zero pressure gradient boundary-layer 
computations with different freestream values for w have 
been performed. It has been shown in [131 that the correct 
freestream values for u outside the boundary-layer are: 


u f = 


4 u 2 


(31) 


where u T = \fr w jp is the friction velocity, defined in terms 
of the wall shear-stress t w , U t is the freestream velocity 
and b* is the displacement thickness. For the first set of 
computations, the above value has been specified at the in- 
flow boundary in the freestream for both the original and 
the BSL k — u> model. Then, the above value was reduced 
by four orders of magnitude and the computations were re- 
peated with both models. Note that the freestream value of k 
was also reduced in order to keep the freestream value of the 
eddy-viscosity constant (equal to the molecular viscosity). 
Figure 3 shows eddy-viscosity profiles for the original and 
the BSL k - w model. The eddy-viscosity of the original 
model changes by almost 100% due to the changes in u f, 
whereas the BSL model gives the same results for both cases. 
The strong sensitivity of the original model to wy is clearly 
unacceptable and can lead to a severe deterioration of the 
results for complex flows, as will be shown later. Results of 
the SST model were also found to be independent of u f. 


If u goes to zero and is finite (typically a fraction of 
the molecular viscosity), the production term for the eddy- 
vicosity goes to infinity for small disturbances in the strain 
rate. A simple way to prevent this from happening is to 
compute both the production term of k, P k , and the dissipa- 
tion term, D k , and than to limit the production term by the 
following formula: 

P k = min(P k , 20- D k ). (30) 

This limiter has been very carefully tested and it was found 
that even for complex flows the ratio of P k j D k reaches max- 
imum levels of only about two inside shear layers. There- 
fore, Equation (30) does not change the solution but only 
eliminates the occurrence of spikes in the eddy-viscosity due 
to numerical "wiggles" in the shear-strain tensor. It also 
eliminated the unphysical buildup of eddy-viscosity in the 
stagnation region of an airfoil, as reported in [14]. Note that 
this is not a specific problem of the k - w model but has also 
been observed with the A: — e model. 


In each of the following comparisons between the dif- 
ferent models, ut for the original k - u model was always 
chosen according to equation (31) whenever possible. In 
cases where the freestream values had to be chosen differ- 
ently, it will be mentioned in the text. Figure 4 shows a 
comparison of the SST, the BSL, the original it - u and the 
Jones-Launder (JL) k - e model (all JL model computations 
have been performed with damping functions as given in 
[19]) for a zero pressure gradient flat plate boundary layer. 
Obviously, all models predict the correct c^-distributions 
and velocity profiles. The k - u models can be run with the 
first gridpoint as far out as y + = 3 without a deterioration 
of the results. 

Free Shear Layers 

For free shear layers the SST and the BSL model reduce 
to the same model (Fj = 0; Fj = 0), which will be called 
new k — u model in this subsection. Note that the new 
k-u model is formally almost identical to the Jones-Launder 
k - e model. However a small cross-diffusion term has 
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been ..ejected in the derivation of equation (4). In order 
to show that this term is truly negligible, the original k — c 
model is also included in the comparison. All computations 
were performed with 200 gridpoints across the layer. A grid 
refinement study on a grid with 300 points gave the same 
results. 

Figure 5 shows the results of the solution of the equi- 
librium far wake equations for the different models. The 
results are compared to the experimental data of Fage and 
Falkner [20], Obviously, the new k - w model and the JL 
k - f model produce almost identical results. The original 
k - w model predicts a somewhat lower spreading-rate than 
the other models. As is well known for two-equation mod- 
els, they fail to predict the smooth behavior of the data at 
the edge of the layer. The freestream value for the original 
k - ^ model has been derived from an expression similar to 
equation (31) [13]. 

Figure 6 shows solutions for a self-similar plane jet, as 
reported in [21], Again, the new k— u> and the JL model pro- 
duce almost identical results and are in very good agreement 
with the experiments, whereas the original k - u> model 
shows a rather peculiar behavior. A major difference be- 
tween the far-wake flow and the present flow is that the 
freestream velocity is zero (still air). The only acceptable 
freestream value for u> in still air is w j - 0, as can be seen 
from the equilibrium equations [13] and from physical intu- 
ition, The specification of small values for w f leads to large 
eddy-viscosities in the original k-ui model, as demonstrated 
in Fig. 3 for the flat plate boundary-layer. The same is true 
in the present flow to an even larger extent, because of the 
missing wall influence. The original k — w model predicts 
about five times as high an eddy-viscosity as the other two 
models, resulting in the large spreading rates shown in Fig. 6. 

A comparison with the free mixing-layer experiment of 
Liepman and Laufer [22] is shown in Fig. 7. Note that the 
freestream velocity below the mixing-layer is zero, leading 
to the same problem with the original it — w model as for the 
plane jet. The other two models again produce almost iden- 
tical results in acceptable agreement with the experiments. 

Adverse Pressure Gradient Flows 

One of the most important aspects of a turbulence model 
for aerodynamic applications is its ability to accurately pre- 
dict adverse pressure gradient boundary-layer flows. It is 
especially important that a model can predict the location of 
flow separation and the displacement effect associated with 
it. The reason is that the viscous-inviscid interaction has 
a strong influence on the overall pressure distribution and 
therefore on the performance of the aerodynamic body. 

The most widely used test case to measure the perfor- 
mance of turbulence models under adverse pressure gradient 
conditions is the flow reported by Samuel and Joubert [23]. 
It is a flat plate boundary-layer, developing under an increas- 
ingly adverse pressure gradient. The upstream Reynolds 


number is 1.7 • 10^m The flow is retarded by the pres- 
sure gradient, but not separated. 

Two different sets of computations have been performed. 
At first, the experimental pressure distribution was specified 
at the outer edge of the computational domain (opposite to 
the wall). Since this is not a very straightforward bound- 
ary condition for a Navier-Stokes method, a second set of 
computations was performed based on the specification of an 
inviscid external streamline. The inviscid streamline y{x) s 
is defined by the preservation of mass: 

y{x)s 

m- J pu eX p.(x,y)dy - const. (32) 
0 

Note that the specification of a streamline does not mean that 
the displacement thickness is prescribed like in an inverse 
boundary-layer method. Both computations produced very 
similar results. The solutions shown here are for the pre- 
scribed streamline, which is thought to be the more consis- 
tent boundary condition from a numerical point of view. The 
eddy-viscosity at the inflow boundary was determined from 
the experimental shear-stress and velocity profiles, the turbu- 
lent kinetic energy k from the requirement k = ( — uV)/ni 
and w and r from the definition of the eddy-viscosity. The 
same grid of 90x90 points as in [12] was used for the com- 
putations. The results are virtually identical to those on a 
60x60 and a 120x120 grid. 

Figure 8 shows a comparison of the computed and the 
experimental skin-friction distribution. All three k - u mod- 
els reproduce the experimental data almost exactly, whereas 
the JL k - f model gives significantly to high values. 

Figure 9 shows the same comparison for the velocity 
profiles. There are no large differences to be found between 
the different models. Only the SST model predicts a some- 
what stronger retardation of the flow near the surface. The 
same behavior has been observed with the Johnson-King 
model in [12], 

Turbulent shear-stress profiles are shown in Fig. 10. 
All models are slightly overpredicting the amount of shear- 
stress, with the SST model closest and the JL model furthest 
away from the data. 

It is obvious that the small differences between the so- 
lutions, especially between the different k - w models, do 
not allow final conclusions about the abilities of the mod- 
els to predict adverse pressure gradient flows. It appears 
that the Samuel-Joubert flow does not pose a sufficiently 
strong challenge to the models to assess their potentials for 
this type of flows. The author has reached a similar con- 
clusion in [12], where a solution of the Johnson-King (JK) 
model has shown, that the model did not significantly depart 
from its equilibrium formulation. It is therefore important to 
test models under more demanding conditions, with stronger 
adverse pressure gradients and possibly separation. The fol- 
lowing flowfield, reported by Driver [15], has proven to be a 
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highly self-consistent and demanding test case, and is there- 
fore strongly recommended for the assessment of turbulence 
models under adverse pressure gradients. 

In Driver’s flow, a turbulent boundary-layer develops 
in the axial direction of a circular cylinder. An adverse 
pressure gradient is imposed by diverging wind tunnel walls 
and suction applied at these walls. The pressure gradient is 
strong enough to cause the flowfield to separate. The inflow 
Reynolds number is 2.8 • 10 5 based on the diameter D of the 
cylinder (140mm). 

The boundary conditions for this flow are similar to 
those used for the Samuel-Joubert flow. Again an inviscid 
streamline is extracted from the experimental velocity pro- 
files and a slip condition is applied along this line. The 
inflow conditions are determined from the experimental pro- 
files in the same way as described above. The computations 
have been performed with a three-dimensional version of 
the code. In order to account for the axial symmetry, three 
closely spaced circumferential planes where introduced and 
symmetry conditions were applied. A 60x3x60 grid [12] 
was used for the present computations. A computation on a 
100x3x100 grid gave almost identical results. 

Figure 1 1 shows the wall pressure distribution for Driver’s 
flow as computed by the different models. The SST model 
gives superior results to the other models due to its ability 
to account for the transport of the principal turbulent shear- 
stress. As expected, the JL k - e model produces the worst 
results, and the BSL and the original k-w model being close 
to each other in the middle. 

Figure 12, depicting the wall shear-stress distribution 
for Driver's flow, shows that the SST model predicts the 
largest amount of separation, whereas the JL model stays 
firmly attached. Again, the BSL and the orig. k - w model 
produce very similar results, in good agreement with the 
experiments. 

The differences between the models can be seen more 
clearly in Fig. 13 which shows the velocity profiles. The 
SST model clearly produces the best agreement with the 
experiments. The larger displacement effect predicted by 
this model is reflected in the flattening of the c p -distribution 
as observed in Fig. 11. The orig. k — w model predicts 
slightly better results than the BSL model, and the JL k - e 
model shows very little sensitivity to the pressure gradient, 
as already reflected in Figs. 1 1 and 12. 

The reasons for the different behavior of the models 
can be seen in the following two pictures. Figure 14 com- 
pares turbulent shear-stress profiles at different stations. The 
experimental data are shown both, in a surface (Carth.) ori- 
ented and in a streamline (Strml.) oriented coordinate sys- 
tem. (Note that the numerical results were, due to the eddy- 
viscosity formulation, not found to be sensitive to small 
changes in direction). 


The JL model obviously predicts significantly higher 
shear-stress levels than the other models, especially in the 
region where separation is approached. This in turn leads to 
the firmly attached velocity profiles of Fig. 13. The differ- 
ences between the models can be seen more clearly by look- 
ing at the eddy-viscosity distributions. Figure 15 shows the 
maximum value of the kinematic eddy-viscosity profiles for 
all x-stations, nondimensionalized by u t 8* . The SST model 
predicts the reduction of this quantity due to the adverse pres- 
sure gradient in very good agreement with the experiments. 
The BSL and the orig. k - w model are very close to each 
other up to separation (around x/D = 0), whereas the orig. 
model is closer to the experiments in the recovery region. 
Both models give consistently to high values for the maxi- 
mum eddy-viscosity in the adverse pressure gradient region. 
The k - e model falls only barely below the value of 0.0168 
recommended by Clauser for equilibrium boundary-layers 
(and used in the Cebeci-Smith model) and thereby fails to 
account for the nonequilibrium effects altogether. 

Backward-Facing Step Flow 

Results for the flow over a backward facing step as 
reported by Driver and Seegmiller [25] will be discussed 
next. This flowfield was a test case in the 1981 Stanford 
conference for the evaluation of turbulence models. How- 
ever, most of the computations at the time were performed on 
comparatively coarse grids and there is substantial evidence 
that significantly finer grids have to be used in order to ob- 
tain grid-independent results [26]. The present computations 
have been performed on a 120x120 grid, with substantial grid 
refinement near the step. Figure 16 shows the distribution 
of gridpoints. As with the other flowfields a grid refinement 
study was made. The present results are virtually identical 
to those performed on a 90x90 and on a 240x240 grid. The 
Reynolds number, based on the upstream momentum thick- 
ness 0 is, Re@ = 5000 and the ratio of the boundary- layer 
thickness to the step height is about 1.5. The expansion ratio 
(height of the tunnel behind the step divided by the height of 
the tunnel in front of the step) is 1.125. 

Figure 17 shows a comparison of computed and ex- 
perimental skin friction distributions. The k - u> models 
all perform significantly better than the k - t model. The 
reattachment length of the four models are 6.4 (SST), 5.9 
(BSL), 6.4 (org. k - w) and 5.5 (JL k - e) compared to 
a value of about 6.4 in the experiments. The reattachment 
length predicted by the k - e model is better than previously 
reported, certainly as a result of the fine grid employed in the 
present computations (see also [26]). However, the model 
predicts significantly too large variations of cy in the recir- 
culation and the reattachment region. The good results of 
the k - u models show that it is not necessary to account for 
the anisotropy of the stress tensor, as suggested in [27], in 
order to get accurate results for the reattachment length. 

The surface pressure distribution, shown in Fig. 18, re- 
flects the trends already seen in c y. The larger the separation 
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region predicted by the model the larger is the displacement 
effect of the boundary layer and the smaller is the pressure 
rise in the expansion region after the step. 

Figure 19 shows a comparison of the velocity profiles. 
All models fail to capture the relaxation downstream of reat- 
tachment correctly. The results of [27] show that this is also 
true for more complex models which account for anisotropy 
effects. 

NACA 4412 Airfoil Flow 

The following set of computations is for the flow around 
a NACA 4412 airfoil at 13.87 degrees angle of attack. The 
Reynolds number with respect to the chord length is Re - 
1.52 10 6 . Experimental data for this flow have been reported 
by Coles and Wadcock [28], The grid for the computations 
consists of 241x61 points and was made available by S. 
Rogers [18], It is similar to the one used in [29], 

In the experiment the transition was fixed by transition 
strips at x/c of 0.023 and 0.1 for the upper and lower sur- 
face respectively. As reported in [29], if these locations are 
specified in the computations, a laminar separation bubble 
appears before the transition point on the upper surface of 
the airfoil. This separation bubble was not observed in the 
experiments which indicates that transition may take place 
already before the strip is reached. Computations have been 
performed with and without a specified transition location 
and differences between the computations are small. Results 
are given here for the case were transition was not specified, 
so that the models picked their own transition location. Tran- 
sition takes place at a downstream station of x /c « 0.006 on 
the upper and x/c ^ 0.06 on the lower surface of the airfoil. 

Figure 20 shows a comparison of the computed and 
the experimental velocity profiles at different streamwise 
stations. The results are similar to those for the separated 
case of Driver, Fig. 13. Again, the SST model predicts the 
displacement effect in very good agreement with the ex- 
periments. The BSL model is showing some response to 
the pressure gradient, and produces results similar to those 
reported in [29] for the Baldwin-Barth model. Another in- 
teresting result of this computation is that the original k - w 
model predicts velocity profiles even further away from the 
experiments than the Jones-Launder k-e model. The reason 
for the poor performance of the orig. k — ui model lies in its 
freestream dependency. In order to understand the problem, 
one has to look at the development of u> from the inflow 
boundary to the leading edge of the airfoil. In this inviscid 
region production and diffusion of 'w are zero. The balance 
in the w equation reduces therefore to: 

U s ^ = -fa 2 (33) 

as 

where s is the streamline direction and U s is the velocity in 
this direction. Assuming that U 3 is constant and equal to 


Uoo, the equation can be solved to give: 

fa*) = —g — p. (34) 

The largest value of w that can be achieved for a certain 
distance s from the inflow boundary is: 

fas) = -J- (35) 

Too s 


corresponding to an infinitely large w 0 c . As s gets large-this 
maximum value becomes small. In the present computations, 
the distance between the inflow boundary and the airfoil is 
about 15 chord lengths. Non-dimensionalizingall quantities 
with Uoo and the chord length, c, leads to a freestream value 
of w in the leading edge region of the airfoil of about ui c = 1 
whereas the formula given in [13] for the estimation of the 
correct freestream value: 


faJUet' 


(36) 


indicates that it should be about three orders of magnitude 
larger. The low freestream value of u> in turn leads to the very 
high eddy-viscosities shown in Fig. 3 which in tum prevent 
separation. This example clearly shows the dangers of using 
the orig. k - u model for aerodynamic applications. If the 
correct freestream values could be specified, the results of 
the orig. k — u should be very close to those of the BSL 
model. 

Figure 21 shows a comparison of the computed and the 
experimental surface pressure distributions. The agreement, 
especially for the SST model is not as good as should be 
expected from the velocity profiles shown in Fig. 20. Al- 
though the SST model predicts the displacement effect of the 
boundary layer almost exactly, it fails to reproduce the pres- 
sure distribution. This indicates an inconsistency between 
the computations and the experiments. Likely candidates 
are blockage effects in the wind tunnel (however, including 
wind tunnel walls in the computations does not improve the 
results [18]) or three-dimensional effects in the experiment. 

Transonic Bump Flow 

The final test case is the axisymmetric transonic shock- 
wave/turbulent boundary layer experiment of Bachalo and 
Johnson [30]. In this experiment, an axisymmetric boundary 
layer interacts with a shock wave created by a circular arc, 
as shown in Fig. 22. It is beyond the scope of this paper 
to present a detailed study of transonic flows and only the 
highest Mach number case will be shown. The Mach num- 
ber for this experiment is 0.925. The number of gridpoints 
used was 150x3x80. Grid-independence was established by 
using different grids (129x3x60 and 180x3x100). Figure 23 
shows the wall pressure distribution computed by the three 
different k — w models, compared with the experiment. The 
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SST model predicts significantly better results than the BSL 
and the original k - u model, due to its improved transport 
features. A compressible version of the k - t model has not 
yet been coded and results for that model are therefore miss- 
ing in the comparison. Earlier results of [31] show however 
a similar behavior as the BSL and the original k - u model. 
Detailed comparisons for transonic flows will be presented 
at a later time. 


Conclusions 

Two new turbulence models have been derived from 
the original k - u model of Wilcox [9]. The motivation 
behind the new baseline (BSL) model was to eliminate the 
freestream dependency of the k —u model but retain its sim- 
ple and reliable form, especially in the near wall region. In 
order to achieve this goal, a switching function was designed 
that can discriminate between the inner part (appr. 50%) of 
a boundary-layer and the rest of the flowfield. In this inner 
part the original k — u model is solved, and in the outer wake 
region a gradual switch to the standard k - e model, in a 
k — u formulation, is performed. 

The BSL model was then used to derivea model that can 
account for the transport of the turbulent shear stress (Shear- 
Stress Transport or SST model). The derivation of the model 
was inspired by the success of the Johnson-King (JK) model. 
The main assumption in the JK model that the principal 
turbulent shear-stress is proportional to the turbulent kinetic 
energy was incorporated into the new SST model. This 
modification ensures that the principal turbulent shear-stress 
satisfies the same transport equation as the turbulent kinetic 
energy. It is designed to act only inside the boundary-layer 
in order to retain the k - e model (transformed to a k—u 
formulation) for free shear-layers. 

Both models were applied to a selection of well docu- 
mented research flows, that are meaningful for aerodynamic 
applications. The results of the computations were compared 
against solutions of the standard k—u and the standard k - e 
model, as well as against experimental data. 

The free shear-layer computations have shown that the 
new models give results almost identical to those of the k — r 
model. Another important aspect of those computations is 
that they show clearly the strong ambiguity in the results of 
the original k-u model with respect to freestream values. 

The central part of the comparisons is for the behavior 
of the models under adverse pressure gradient conditions. 
The computations of the Samuel-Joubert flow, as well as 
Driver’s separated adverse pressure gradient flow show that 
the SST model gives highly accurate results for this type of 
problem. The BSL and the original k - u model produce 
rather similar results, provided the correct freestream values 
are specified in the latter. 

Computations were also performed for the backward 
facing step flow of Driver and Seegmiller [25]. A very fine 


grid was employed to ensure grid independence of the re- 
sults. For this problem, the original k - u and the SST 
model give very accurate results. They predict the reattach- 
ment length within the uncertainty of the measurements and 
give an accurate representation of the wall pressure distribu- 
tion. The BSL model gives about 8% too small values for 
the reattachment length. These results are still very accurate, 
considering the notorious difficulties this flow poses to nu- 
merical assessment. All models fail to predict the relaxation 
of the velocity profiles downstream of the reattachment point 
correctly. 

Computations for a NACA 4412 airfoil at an angle of 
attack near maximum lift condition confirm the findings of 
the adverse pressure gradient computations. The SST model 
predicts highly accurate velocity profiles, almost identical 
to those of the experiments. The BSL model has a smaller 
sensitivity to the adverse pressure gradient than the SST 
model and therefore predicts less retarded profiles. A very 
surprising result of the computations is that the original k-u 
model gives even less accurate solutions than the Jones- 
Launder k — < model. The reason for the failure of the 
model is again its freestream dependency. This computation 
clearly shows that the original k -u model cannot be applied 
unambiguously for industrial applications. 

The last set of computations is for a transonic shock- 
wave/turbulent boundary layer interaction. The accurate 
prediction of the shock location by the SST model shows 
that the good performance of this model for incompressible 
applications can be extended to transonic flows. 

Appendix 

The Baseline (BSL) Model 


Dpk du, , 5 


dk 


(A- 1) 


Dpu 7 du, 2 . d 
Dt v t T ' ] dxj 0pU + dx } 

, on r . 1 dk du 

+ 2(1 - Fl )pa u2 - — — 


{p + <JuHt) 


1 . 
(A-2) 


The constants <p of the new model are calculated from 
the constants, 4 >\ , <£ 2 . as follows: 


4> = F^\ +(1 - Fi)<£> 2 - (A-3) 

The constants of Set 1 (^j) are (Wilcox): 

°k\ = 0.5, <r wl = 0.5, /?! = 0.0750, (A-4) 

/?* = 0.09, k = 0.41, 71 = /?!//?* - <7^ k 2 /^. 
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The constants of Set 2 (f 2 ) are (standard k - e): 


and the eddy-viscosity is defined as: 


a k2 = 1.0, cr w2 = 0.856, 0 2 = 0.0828, (A-5) 

p* = 0.09. k = 0.41 , 72 = &2 /&* ~ Vul* 2 /'/$*■ 


The function F j is defined as follows: 


F\ = tanh(arg^) (A-6) 


with: 


" t 


a 1 k 

max(ai<j; flF 2 ) 


(A- 14) 


where ft is the absolute value of the vorticity. F 2 is given 
by: 

F 2 — tanh(arg^) (A- 15) 


with: 


arg 2 — max 


(2 


\/k 

0.09u V 


500 u 

2 ’ 
y ^ 


(A- 16) 


. ( . Vk 

arg { = nnn(mat(— - 


500i/ Aprr^k \ 
CD kujy" 


(A-7) 


where y is the distance to the next surface and C D kuJ is the 
cross-diffusion term of equation (A-2): 


Important detail !: 

In applying this model, it is important that the reader is 
aware of the following ambiguity in the formulation of the 
production term of ui for the SST model. The definition of 
the production term of w is sometimes written as: 


1 dk dt 


CD kul = max(2pcr 2 — -\ 

u axj oxj 




(A-8) 


Plj = 7 


’.jj du, 

k IJ dxj 


(A- 17) 


The eddy-viscosity is defined as: 

v t = - (A-9) 

JJ 

and the turbulent stress tensor r tJ = -pu'u' is: 


which introduces the nondimensional group vtj > n front of 
the strain rate tensor. In the original and in the BSL model 
this group is equal to one and the two formulations for P L 
are therefore identical. This is not the case for the SST 
model because of equation (A- 14). The SST model has been 
calibrated with respect to equation (A-2) and equation (A- 1 7) 
should therefore not be used. 


^ = + < A -«» 


dx l 3 ctej. t} ‘ 


The following choice of freestream values is recommended: 


Woe = (1 — ► 10)— p- (A-ll) 

r'foo = 10 -3 £/ 

kec = r'too^oo (A- 12) 


where L is the approximate length of the computational do- 
main. 


Appendix 

The Shear-Stress Transport (SST) Model 

The SST model is identical to the above formulation, 
except that the constants, <j>\, have to be changed to: 
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(a) Function F j (b) Function Fi 

Fig. 1 Blending functions F\ and FS versus y/S for different velocity profiles. 
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Fig. 3 Freestream-dependency of the eddy-viscosity for the original and the BSL k —u model 
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(a) Skin friction 


(b) Law-of-the-wall 


Fig. 4 Flat plate boundary-layer. 
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Fig. 5 Comparison of the velocity distribution of three different two-equation models for the far wake flow of Fage and 
Falkner 
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Fig. 6 Comparison of the velocity distribution of three different two-equation models for the plane jet flow of Wygnanski 
and Fiedler 
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Fig. 7 Comparison of the velocity distribution of three different two-equation models for the free mixing-layer of Liepman 
and Laufer 



Fig. 8 Wall shear-stress distribution for Samuel-Joubert flow. 



Fig. 9 Velocity profiles for Samuel-Joubert flow at x= 1.16, 1.76, 2.26, 2.87, 3.40 m. 
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Fig. 14 Turbulent shear-stress profiles for Driver’s adverse pressure-gradient flow at x/D=-0.091, 0.363, 1.088, 1.633, 





















Fig. 21 Surface pressure distribution for a NACA 4412 airfoil at 13.87 degrees angle of attack. 
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Fig. 22 Transonic bump flow experiment of Bachalo and Johnson 




Fig. 23 Comparison of cp-distributions for transonic bump flow at M=0.925 





